Preparations

Load the necessary libraries

library(car)
library(rstanarm) # for fitting models in STAN
library(brms) # for fitting models in STAN
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(DHARMa) # for residual diagnostics
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(tidyverse) # for data wrangling etc
library(patchwork) # for multiple plots
library(ggridges) # for ridge plots

Scenario

Loyn (1987) modelled the abundance of forest birds with six predictor variables (patch area, distance to nearest patch, distance to nearest larger patch, grazing intensity, altitude and years since the patch had been isolated).

Regent honeyeater

Format of loyn.csv data file

ABUND DIST LDIST AREA GRAZE ALT YR.ISOL
.. .. .. .. .. .. ..
ABUND Abundance of forest birds in patch- response variable
DIST Distance to nearest patch - predictor variable
LDIST Distance to nearest larger patch - predictor variable
AREA Size of the patch - predictor variable
GRAZE Grazing intensity (1 to 5, representing light to heavy) - predictor variable
ALT Altitude - predictor variable
YR.ISOL Number of years since the patch was isolated - predictor variable

The aim of the analysis is to investigate the effects of a range of predictors on the abundance of forest birds.

Read in the data

loyn <- read_csv("../public/data/loyn.csv", trim_ws = TRUE)
glimpse(loyn)
## Rows: 56
## Columns: 7
## $ ABUND   <dbl> 5.3, 2.0, 1.5, 17.1, 13.8, 14.1, 3.8, 2.2, 3.3, 3.0, 27.6, 1.8…
## $ AREA    <dbl> 0.1, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.…
## $ YR.ISOL <dbl> 1968, 1920, 1900, 1966, 1918, 1965, 1955, 1920, 1965, 1900, 19…
## $ DIST    <dbl> 39, 234, 104, 66, 246, 234, 467, 284, 156, 311, 66, 93, 39, 40…
## $ LDIST   <dbl> 39, 234, 311, 66, 246, 285, 467, 1829, 156, 571, 332, 93, 39, …
## $ GRAZE   <dbl> 2, 5, 5, 3, 5, 3, 5, 5, 4, 5, 3, 5, 2, 1, 5, 5, 3, 3, 3, 2, 2,…
## $ ALT     <dbl> 160, 60, 140, 160, 140, 130, 90, 60, 130, 130, 210, 160, 210, …

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ log(\mu_i) = \boldsymbol{\beta} \bf{X_i} \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the additive effects of the scaled versions of distance (ln), distance to the nearest large patch (ln), patch area (ln), grazing intensity, year of isolation and altitude on the abundance of forest birds.

Fit the model

MCMC sampling diagnostics

Model validation

Partial effects plots

Model investigation

Further analyses

Contrasts

brms

loyn.list <- with(loyn, list(AREA = c(min(AREA), mean(AREA), max(AREA))))

newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  as.data.frame()
head(newdata)
loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list) %>%
  regrid() %>%
  pairs()
## AREA =    0.1:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2       9.878     0.635     19.62
##  1 - 3      10.664     1.222     20.24
##  1 - 4      16.464     6.831     25.25
##  1 - 5      17.421     9.546     25.95
##  2 - 3       0.847    -6.897      8.37
##  2 - 4       6.819    -1.114     13.95
##  2 - 5       7.541     1.409     14.27
##  3 - 4       5.917    -2.054     12.85
##  3 - 5       6.720     1.353     13.25
##  4 - 5       0.627    -4.391      7.33
## 
## AREA =   69.3:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2      -3.173   -10.984      4.27
##  1 - 3      -5.151   -15.230      2.86
##  1 - 4     -25.969   -66.460      3.54
##  1 - 5       8.640   -12.511     23.70
##  2 - 3      -2.012   -12.553      8.75
##  2 - 4     -22.585   -67.948      5.13
##  2 - 5      11.689    -9.172     29.49
##  3 - 4     -20.470   -63.429      8.65
##  3 - 5      14.107    -9.102     31.25
##  4 - 5      33.848    -1.088     81.45
## 
## AREA = 1771.0:
##  contrast estimate lower.HPD upper.HPD
##  1 - 2     -22.681   -58.074      4.92
##  1 - 3     -30.888   -79.271      4.95
##  1 - 4    -214.126 -1096.279     16.27
##  1 - 5     -29.823  -326.031     37.12
##  2 - 3      -8.558   -66.073     40.75
##  2 - 4    -190.594 -1092.333     46.21
##  2 - 5      -6.397  -321.484     82.77
##  3 - 4    -181.636 -1072.676     67.70
##  3 - 5       2.164  -310.594    111.55
##  4 - 5     164.496  -504.042   1206.70
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95
newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  gather_emmeans_draws()

newdata %>%
  median_hdci() %>%
  ggplot() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(y = .value, ymin = .lower, ymax = .upper, x = contrast)) +
  facet_wrap(~AREA) +
  coord_flip()

emmeans(loyn.brm4b, ~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  gather_emmeans_draws()
newdata.p <- newdata %>% summarise(P = sum(.value > 1) / n())
g <- newdata %>%
  ggplot() +
  geom_vline(xintercept = 1, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = contrast,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
  facet_grid(~ round(AREA, 1))

# g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = round(P,3)))
# g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = paste('P = ',round(P,3))), hjust = -0.2, position = position_nudge(y = 0.5))

Summary figure

References

Loyn, R. H. 1987. “Nature Conservation: The Role of Remnants of Native Vegetation.” In, edited by D. A. Saunders, G. W. Arnold, A. A. Burbridge, and A. J. M. Hopkins. Chipping Norton, NSW: Surrey Beatty & Sons.